Giacomo Moretti

Luca Zaccarian

Franco Blanchini

Quadratic constrained periodic optimisation for bandlimited linear systems via the Fourier-based method

Keywords: Fourier series, frequency domain, dynamic programming, convex optimisation, periodic, bandlimited

come

Quadratic constrained periodic optimisation for bandlimited linear systems via the Fourier-based method

Giacomo Moretti, Luca Zaccarian, Franco

Introduction

Dynamic programming techniques aim at determining optimal inputs for a dynamical systems trajectory tracking while minimising a performance index (namely, a cost) in the presence of constraints [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF]. Since analytical solutions are limited to a small set of practical cases, optimal control problems are most often tackled numerically.

In this regard, two main classes of numerical methods for optimal control exist [START_REF] Rao | A survey of numerical methods for optimal control[END_REF], namely, time-domain methods (also called time-marching methods) and functional-approximation methods (also called collocation methods). Time-domain methods rely on the approximation of a continuous-time plant as a discrete system. According to these approaches, the continuous-time control input is approximated by a set of evaluations at a finite set of time instants. These evaluations are hence used as the unknown of a finite-dimensional minimisation problem. Functional-approximation methods, in contrast, rely on the representation of the system state/input as a linear combination of independent functions (either defined globally on the optimisation horizon, or locally on a set of time intervals), enforcing the constraints at a set of intermediate points and reducing the optimisation problem to the calculation of the combination coefficients [START_REF] Bellman | Functional approximations and dynamic programming[END_REF]. As compared to time domain approaches, functional-approximation methods allow representing continuous-time control functions through a reduced number of parameters, hence limiting the computational burden and possibly rejecting non-smooth control solutions, typical of time domain approaches.

A relevant class of dynamic programming problems concerns optimal control of steady-state oscillating systems in the periodic regime. These problems are of interest in different engineering fields, such as energy harvesting from oscillating sources (human walking, structural vibrations, ocean wave energy) via mechanical oscillators [START_REF] Mitcheson | Energy harvesting from human and machine motion for wireless electronic devices[END_REF], and driving of electrical circuits operating in alternated current [START_REF] Torquati | A hybrid predictive control algorithm for tracking in a single-phase DC/AC inverter[END_REF]. This paper deals with a theoretical framework for dynamic programming problems in steady-state periodic 1 systems via a so-called Fourier-based method [START_REF] Nagurka | Fourier-based optimal control of nonlinear dynamic systems[END_REF], i.e. a functional-approximation approach using a finite Fourier series (FS) expansion to represent the system states and inputs [START_REF] Hussaini | Spectral collocation methods[END_REF]. This method matches particularly well optimal control problems in dynamic periodic systems, as it provides an intrinsically periodic representation of the states and of the inputs with no need for additional periodicity constraints. The application of functional spectral methods based on the FS for optimisation problems is limited to a few works. Yen & Nagurka [START_REF] Yen | Fourier-based state parameterization for linear quadratic optimal control[END_REF][START_REF] Nagurka | Fourier-based optimal control of nonlinear dynamic systems[END_REF][START_REF] Yen | Linear quadratic optimal control via Fourier-based state parameterization[END_REF] proposed for the first time a framework for optimisation problems with Fourier-based state parametrisation, first with reference to linear quadratic (LQ) problems [START_REF] Yen | Fourier-based state parameterization for linear quadratic optimal control[END_REF][START_REF] Yen | Linear quadratic optimal control via Fourier-based state parameterization[END_REF], and then to the general case of systems with nonlinear dynamics [START_REF] Nagurka | Fourier-based optimal control of nonlinear dynamic systems[END_REF]. Recently, Bacelli et al. [START_REF] Bacelli | A control system for a self-reacting point absorber wave energy converter subject to constraints[END_REF][START_REF] Bacelli | Numerical optimal control of wave energy converters[END_REF] applied Fourier-based transcription for the identification of optimal control strategies for ocean wave energy converters, aimed at maximising the harvested energy in the presence of operational constraints. In those works [START_REF] Yen | Linear quadratic optimal control via Fourier-based state parameterization[END_REF][START_REF] Bacelli | Numerical optimal control of wave energy converters[END_REF], continuous-time constraints are reformulated in a relaxed manner and enforced only at a finite number of points, generally leading to solutions which are not strictly feasible for the original problem.

In this paper, we propose a Fourier-based framework that specifically applies to the case of periodic oscillating systems. As in [START_REF] Yen | Fourier-based state parameterization for linear quadratic optimal control[END_REF][START_REF] Yen | Linear quadratic optimal control via Fourier-based state parameterization[END_REF], we address LQ problems involving the minimisation of a quadratic cost in a system subject to hard constraints with linear dynamics, leading to nonlinear optimal responses. We first show that some optimisation problems of technical interest in different fields (namely, energy harvesting from vibrating sources and control of inverters) can be reduced, in a first approximation and for the sake of landscaping analysis, to a same LQ periodic convex formulation. With the aim of realistically describing technical applications, we then set the attention on systems with smooth bandlimited excitation, input and response. In this regard, we show that the optimisation problem can be cast using a finite-dimensional FS description of the system, leading to a convex quadratic formulation. Finally, we propose a new algorithm for the identification of a sub-optimal strictly feasible solution (respecting the constraint over the entire continuous time horizon) leading to an optimal cost arbitrarily close to the optimal cost for the original problem.

The paper is organised as follows. Sect. 2 introduces the class of LQ problems under investigation and shows that a diversity of practically relevant problems belong to that class. Under the assumption of bandlimited periodic input and excitation, Sect. 3 presents a compact transcription of the system dynamics, the cost function and the constraints in terms of truncated FS, leading to a reformulation of the original LQ problem as a convex quadratic minimisation problem with parametric constraints.

Sect. 4 provides a practical algorithm for the identification of a nearly-optimal feasible solution for the obtained minimisation problem through a finite-dimensional discretisation of the constraints. Sect. 5 illustrates the application of the proposed framework on a few numerical case studies. Finally, Sect. 6 draws the conclusions.

2 Linear quadratic optimisation of steady-state periodic responses

We consider a class of LQ optimal control problems [START_REF] Han | A unified numerical scheme for linearquadratic optimal control problems with joint control and state constraints[END_REF] involving a steady-state oscillating linear dynamical system subject to a periodic excitation, a cost represented by a quadratic function of the system inputs and states, and linear inequality constraints on the state and/or the input. The dynamical system is asymptotically stable and described by:

ẋ(t) = Ax(t) + B(u(t) + d(t)), (1)
where t ≥ 0 is the continuous time, x ∈ R m is a state vector, u ∈ R p is a vector of p control inputs, d ∈ R p is a known input-matched bounded periodic excitation with period T , namely, d(t) = d(t + T ) for all t ≥ 0, A ∈ R m×m is a Hurwitz matrix called the system matrix, and B ∈ R m×p is called the input matrix. We consider LQ optimisation problems in the variable u involving a steady-state oscillatory response of the system. We thus restrict the analysis to a time interval with length T and functions d, u and x belonging to the space L 2 ([0, T ]), namely square integrable functions on [0, T ]. The optimisation problem that we address is defined below.

Problem 1 Continuous-time LQ periodic optimisation problem:

           min u∈L 2 ([0,T ]) J(u) = T 0 (x(t) P x(t) + u(t) Qu(t)) dt subject to (1) A x x(t) + A u u(t) + g(t) ≤ b, ∀t ∈ [0, T ] x(0) = x(T ), (2)
where P = P ∈ R m×m and Q = Q ∈ R p×p are symmetric positive-semidefinite matrices; A x ∈ R l×m and A u ∈ R l×p , g : R → R l is a known, bounded and T -periodic function of time, and b ∈ R l >0 is a vector with positive entries. Parameters A x , A u , g and b define a set of l linear inequalities that bound a feasible set for the problem. It is hereby assumed that such a feasible set has non-empty interior. Since the constraint inequality has two intercepts (g(t) and b), b can be chosen to be positive-valued with no loss of generality. This choice simplifies the use of restrictions of the feasible set that will be exploited later on in the paper. Problem 1 includes a terminal state constraint requiring the periodicity of x. Periodicity of x and d automatically guarantees periodicity of u. Since x and u are related through the linear differential relationship (1), the integrands of both costs can be regarded as quadratic functions of u. Since steadystate oscillations are investigated, no final-state cost is here considered. The cost J and the constraints are convex, therefore the considered problem is a convex optimisation problem.

The considered problem matches a number of physical cases of practical interest, which can be reshaped into the general convex form of (2) in Problem 1.

In the following, we discuss two applications that fit within the framework of Problem 1. The first application regards energy harvesting from a mechanical oscillator, while the second one regards the control of a single-phase inverter.

Application 1. Maximal energy harvesting from a harmonic oscillator with periodic excitation

The first problem under investigation concerns energy harvesting from pulsating sources of mechanical energy. Mechanical energy harvesting is addressed in a considerable number of recent research works, focused on different target applications at different scales and operating frequencies. Applications include energy harvesting from vibrations (generated, e.g., by industrial machines) for self-powered wireless sensors [START_REF] Beeby | Energy harvesting vibration sources for microsystems applications[END_REF]; energy harvesting from human walk to recharge low-power portable electronics [START_REF] Mitcheson | Energy harvesting from human and machine motion for wireless electronic devices[END_REF]; and scavenging of sea waves power for large-scale electricity production [START_REF] Falcao | Wave energy utilization: A review of the technologies[END_REF].

System model

An oscillating mechanical energy harvester can be modelled as a mass-spring-damper system subject to a periodic excitation, as schematically shown in Fig. 1a. The excitation is provided either by the periodic displacement of an inertial frame moved by external vibrations (Fig. 1a left), or by an excitation force directly applied on the mass (Fig. 1a right). Energy extraction is performed by a controllable power take-off (PTO) transducer. The first case (Fig. 1a left), referred to as inertial energy harvester, suitably models energy scavenging devices from vibrations and human motion [START_REF] Beeby | Energy harvesting vibration sources for microsystems applications[END_REF], in which the displacement of the frame holding the device (e.g., a vibrating structure, a moving human body surface) is assumed to be known and is not substantially altered by the dynamics of the seismic mass. The second model (Fig. 1a right) roughly represents the operation of a sea wave energy converter (WEC), consisting in a floating oscillating body subject to wave-induced loads [START_REF] Moretti | Resonant wave energy harvester based on dielectric elastomer generator[END_REF]. In practice, the real dynamics of a WEC include a more sophisticated set of interactions (the system in Fig. 1a-right only represents an approximation), including frequency-dependent dynamical parameters and nonlinear hydrodynamic contributions [START_REF] Li | A synthesis of numerical methods for modeling wave energy converter-point absorbers[END_REF].

With reference to the systems shown in Fig. 1a, we call M the inertia of the oscillator, k the stiffness associated to the restoring forces (e.g., compliance of the joint between the oscillator and the frame in inertial generators, hydrostatic forces in WECs), c the linearised external damping due to dissipative forces. Coefficients M, c and k are positive. We indicate with u the controllable PTO load. Indicating with z the displacement of the mass from its rest position, the following dynamic equation applies:

M z(t) + c ż(t) + kz(t) = u(t) + e(t), (3)
where e is an external excitation. In inertial generators (Fig. 1a left), e depends on the displacement y of the vibrating frame: e = c ẏ(t) + ky(t). In direct-force generators (Fig. 1a right), e simply represents the external excitation force. We hereby consider periodic time-series for the excitation e and the steady-state response of the system. This is a natural assumption in applications involving inertial oscillators: periodic displacements y(t) of the moving frame (associated, e.g., with structural vibrations or human walk) indeed lead to a periodic profile for e. In the case of sea WECs, in contrast, the excitation force due to real waves is not periodic and can be regarded as a superposition of harmonic contributions with random phases. Nonetheless, the analysis of the WEC response and performance in regular sinusoidal waves is a widely used tool for preliminary evaluations, and optimal control in regular waves is used to extrapolate control heuristics for application in real conditions [START_REF] Rosati Papini | Control of an oscillating water column wave energy converter based on dielectric elastomer generator[END_REF].

Optimisation problem formulation

A relevant optimisation goal is to maximise the harvested energy in the presence of a given periodic excitation. Indicating with Ẇu = -u ż the instantaneous power output of the PTO ( Ẇu > 0 if the PTO is positively harvesting energy from the oscillating system), we consider the problem of optimising the total harvested energy over one period in steady-state conditions, in the presence of bounds on the maximum admissible oscillation amplitude z of the moving mass and ū of the PTO load excitation:

               min u∈L 2 ([0,T ]) J(u) = T 0 żudt subject to (3) |u(t)| ≤ ū, ∀t ∈ [0, T ] |z(t)| ≤ z, ∀t ∈ [0, T ] z(0) = z(T ) (4)
Problem (4) has been largely investigated in the field of energy harvesting [START_REF] Renno | On the optimal energy harvesting from a vibration source[END_REF][START_REF] Falnes | Ocean waves and oscillating systems: linear interactions including wave-energy extraction[END_REF]. In the case of monochromatic harmonic excitation, if no constraints ū on the PTO load and z on the system oscillation amplitude are present, energy transfer from the input mechanical energy source to the PTO is maximized when the system response is resonant, i.e., when the system velocity ż is in phase with the excitation d [START_REF] Renno | On the optimal energy harvesting from a vibration source[END_REF].

The presence of constraints leads to nonlinear dynamics of the harvester with latched oscillatory behaviours [START_REF] Halvorsen | Architecture-independent power bound for vibration energy harvesters[END_REF], whose solution has been tackled, in the past, using time domain dynamic programming approaches [START_REF] Hals | Constrained optimal control of a heaving buoy wave-energy converter[END_REF] and Fourier-based approaches [START_REF] Bacelli | Numerical optimal control of wave energy converters[END_REF]. The instantaneous power Ẇu associated with the optimal control u * is not always positive, i.e., the PTO is required to supply power to the oscillator during certain time intervals, leading to some amount of reactive power to be injected in the system through the PTO. In general, the optimal solution u * to (4) is non-smooth.

Problem (4) has an objective function which is, in general, non convex. A convex reformulation of (4) can be obtained in the frequency domain for the case of harmonic excitation e [START_REF] Renno | On the optimal energy harvesting from a vibration source[END_REF][START_REF] Falnes | Ocean waves and oscillating systems: linear interactions including wave-energy extraction[END_REF]. In the following, we consider the case of a generic periodic excitation and we show that the problem can be conveniently resorted to a convex formulation. In particular, it can be cast in the form of (2) in Problem 1 through an appropriate selection of the state x.

Proposition 1 Problem (4) can be cast in terms of a convex quadratic cost.

We consider scalar equation (3) for the system dynamics. Multiplying both members of the equation by ż and computing the integral over one time period allows rewriting the objective function of (4) as follows:

J(u) = T 0 żudt = T 0 c ż2 -że dt = T 0 c ż - e 2c 2 - e 2 4 dt, (5)
where the integrals of the terms M z ż and kz ż have been omitted since they are clearly equal to zero due to the periodicity of z(t). Being e a fixed bounded function of time, minimisation of J in ( 5) is equivalent to the minimisation of

J(u) = T 0 ż - e 2c 2 dt, (6)
which is a convex function of ż. Due to the linear differential relationship between z and u, cost J is also as a convex function of u, hence proving Proposition 1.

Eq. ( 6) clearly shows that, in the absence of constraints, the optimal velocity profile that maximises energy harvesting is given by ż * = e/(2c). In this condition, also known as resonance condition, the mass velocity is in phase with the excitation force.

Proposition 2 Problem (4) can be expressed in the general form (2) of Problem 1.

To prove Proposition 2 we introduce a convenient definition for the system state as follows:

x = ż -e/(2c) z -e p /(2c) , (7)
where e p is any primitive of e (namely, ėp = e). In practice, the state of the system is defined here as the difference between the physical variables (z, ż) and their optimal unconstrained trends. With choice [START_REF] Hussaini | Spectral collocation methods[END_REF] for the state x, Eq. ( 3) can be reduced to a differential equation of the form (1) with m = 2 and a one-dimensional input (p = 1), with the following selections for the system and input matrices:

A = -c/M -k/M 1 0 , B = 1 0 , (8)
where u represents the PTO force and the following expression is used for the external excitation:

d = e 2 - M ė 2c - ke p 2c . (9)
Since e(t) is periodic, d is also periodic.

Based on [START_REF] Nagurka | Fourier-based optimal control of nonlinear dynamic systems[END_REF], optimisation problem (4) can be thus rewritten in the following LQ form:

                   min u∈L 2 ([0,T ]) J(u) = T 0
x 2 1 dt subject to (1), ( 8), ( 9)

|u(t)| ≤ ū, ∀t ∈ [0, T ] x 2 + e p (t) 2c ≤ z, ∀t ∈ [0, T ] x(0) = x(T ) (10)
Problem [START_REF] Bacelli | A control system for a self-reacting point absorber wave energy converter subject to constraints[END_REF] corresponds to Problem 1 (hence proving Proposition 2) with the following equalities:

P = 1 0 0 0 , Q = 0, A x = 0 0 0 0 0 0 1 -1 , A u = 1 -1 0 0 , b = ū ū z z , g(t) = 0 0 e p (t) 2c - e p (t) 2c . (11)

Application 2. Optimal tracking of a reference output voltage in DC/AC converters

The second problem regards optimal driving of DC/AC power converters, namely inverters. Although the technology and control of inverters dates back to decades ago, a quest for improved regulation 140 strategies is currently taking place as a result of new power resources intrinsically operating in DC (such as photovoltaic systems) being connected to the grid [START_REF] Prodanovic | Control and filter design of three-phase inverters for high power quality grid connection[END_REF]. Grid-tied inverters are required to comply with stringent requirements in terms of power output quality. For this reason, the identification of robust regulation strategies, providing limited harmonic distortion, is currently a relevant research topic [START_REF] Torquati | A hybrid predictive control algorithm for tracking in a single-phase DC/AC inverter[END_REF][START_REF] Mariéthoz | Explicit model-predictive control of a PWM inverter with an LCL filter[END_REF].

System model

We hereby consider a single-phase inverter, as schematically shown in Fig. 1b. The device consists in a linear RLC filter in cascade with a bridge of switching components (namely, an H-bridge). We call R, L and C the positive values of the resistance, inductance and capacitance in the RLC branch, respectively. A DC input voltage withconstant amplitude ū is fed into the device on the H-bridge side and converted into an average voltage profile u (with |u(t)| ≤ ū ∀t ∈ [0, T ]), supplied to the RLC through a pulse width modulation (PWM) switching logic. The voltage drop over the capacitor is the inverter output, denoted v. In practical applications, v is required to track a target signal (e.g., a sinusoidal voltage with a frequency equal to the grid frequency). We set our attention on the dynamics of the RLC filter, treating the voltage u as the input and v as the output. The following dynamical model described the circuit of Fig. 1b:

Lq(t) + R q(t) + q(t)/C = u, i(t) = q(t), v(t) = q(t)/C, ( 12
)
where q is the charge on the capacitor and i is the current in the loop.

Optimisation problem formulation

The optimisation problem that we address for [START_REF] Han | A unified numerical scheme for linearquadratic optimal control problems with joint control and state constraints[END_REF] regards the minimization of the difference between the output voltage v(t) and a target periodic trajectory v r (t). Since, in practice, the inverter is controlled acting on the H-bridge switches, the problem can be tackled using two alternative approaches, as observed in [START_REF] Mariéthoz | Explicit model-predictive control of a PWM inverter with an LCL filter[END_REF]. An approach consists in employing the discrete states of the switches as the problem variables. This monolithic approach provides a direct solution for the optimal sequence of the physical states of the switches, but has rather high computational burden since it involves a finite-state formulation that should also incorporate technological constraints on dwell time between consecutive switches. The second approach consists in initially restricting the attention to the average current behaviour, using the profile of an electrical variable (e.g., the RLC filter input voltage) as the unknown of the optimisation problem. The solution to this problem can be employed, in a second phase, as a reference evolution to identify an appropriate switching strategy, e.g. by resorting to a hybrid formulation for the circuit dynamics [START_REF] Torquati | A hybrid predictive control algorithm for tracking in a single-phase DC/AC inverter[END_REF].

We hereby consider the first step of the second approach, and use the input voltage u of the RLC filter as the unknown of the problem. The magnitude of this voltage is clearly required to be smaller than the DC input voltage. The resulting optimisation problem is given by ( 13) below, where the cost function is defined as the squared L 2 norm of the difference between v and v r .

           min u∈L 2 ([0,T ]) J(u) = T 0 [v(t) -v r (t)] 2 dt subject to (12) |u(t)| ≤ ū, ∀t ∈ [0, T ] v(0) = v(T ) (13)
The periodicity of v automatically guarantees that i and v are also periodic.

It is worth noticing that a perfect tracking of the reference trajectory (namely, v * = v r ) is obtained as the solution of (13) whenever the following condition is satisfied:

|Lv r + R vr + Cv r | ≤ C ū (14)
Condition [START_REF] Falcao | Wave energy utilization: A review of the technologies[END_REF] guarantees that the input amplitude constraints are not active. Then the optimal solution coincides with that of the unconstrained scenario. Similarly to Application 1 in Sect. 2.1, in the presence of constraints, the solution to (13) generally leads to non-smooth dynamics.

The following result holds for the general case, regardless of ( 14).

Proposition 3 Problem (13) can be expressed in the LQ form given by Problem 1.

We hereby define a set of state variables representing the tracking error of the state variables with respect to a reference dynamics, in the same fashion as in Sect. 2.1.2. Similarly to [START_REF] Torquati | A hybrid predictive control algorithm for tracking in a single-phase DC/AC inverter[END_REF], in addition to the reference voltage v r on the capacitor, we define a corresponding reference current i r = C vr and we define the state of the system as x = [i -i r , v -v r ] T . Rewriting [START_REF] Han | A unified numerical scheme for linearquadratic optimal control problems with joint control and state constraints[END_REF] in terms of x leads to a linear model of the form (1), with the following definitions for the system matrix, the input matrix and the excitation signal:

A = -R -1 1/C 0 , B = 1/L 0 , d = (LCω 2 0 -1)v r -Ri r . (15)
Since v r and i r are periodic functions, d is also periodic.

155

Using the above-defined error x, optimisation problem (13) takes the following form:

           min u∈L 2 ([0,T ]) J(u) = T 0 x 2 2 dt subject to (1), (15) |u(t)| ≤ ū, ∀t ∈ [0, T ] x(0) = x(T ). ( 16
)
The problem, in this form, clearly belongs to the class of problems defined by (2) in Problem 1 (hence proving Proposition 3), with the following definitions:

P = 0 0 0 1 , Q = 0, A x = 0 2×2 , A u = 1 -1 , g(t) = 0 2×1 , b = ū ū , ( 17
)
where 0 i×j indicates a matrix of zeros in R i×j .

3 Fourier-based dynamic programming for bandlimited periodic systems

The solution to (2) in Problem 1 generally leads to non-smooth periodic dynamics, such that the optimal input and the state evolution include high-frequency harmonics, as it can be verified by applying Pontryagin's maximum principle. Practical applications typically involve smooth bandlimited excitations d and control inputs u, compatible with the physical mechanisms or the machinery through which they are generated. Owing to linear model [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF], this results in a smooth bandlimited evolution of the state x. Instead of pursuing a Pontryagin-based approach or a time-domain numerical approach, which would provide non-smooth solutions, in this section we discuss a methodology to cast the optimisation problems discussed in the previous section in a bandlimited fashion, using a FS approach.

The idea behind the proposed method is to express u as a FS with a finite number of harmonics and to recast the original optimisation problem as a finite-dimensional minimisation problem holding the Fourier coefficients of u as the unknown variables. With reference to linear dynamics (1), we call u n , d n and x n the parametrisation of u, d and x as a set of FS with a finite number n of harmonics of the fundamental frequency ω 0 = 2π/T , where T is the period in optimisation (2):

u(t) ≈ u n (t) = c u 0 + n k=1 c u k cos (kω 0 t) + s u k sin (kω 0 t) , d(t) ≈ d n (t) = c d 0 + n k=1 c d k cos (kω 0 t) + s d k sin (kω 0 t) , x(t) ≈ x n (t) = c x 0 + n k=1 c x k cos (kω 0 t) + s x k sin (kω 0 t) , (18)
where

c u k , s u k , c d k , s d k ∈ R p , c x k , s x k ∈ R m are Fourier coefficients.
Substituting the previous expansions in [START_REF] Bertsekas | Dynamic programming and optimal control[END_REF] provides

Ac x 0 + B(c u 0 + c d 0 ) + n k=1 Ac x k -kω 0 s x k + B(c u k + c d k ) cos (kω 0 t) + As x k + kω 0 c x k + B(s u k + s d k ) sin (kω 0 t) = 0 m×1 . (19)
Due to the orthogonality of the Fourier basis functions, for [START_REF] Falnes | Ocean waves and oscillating systems: linear interactions including wave-energy extraction[END_REF] to hold for all t ∈ [0, T ], each one of the terms in square brackets must be equal to 0 m×1 , leading to the following set of linear equalities (which must hold for all k = 1, . . . , n):

-A M0 c x 0 = B(c u 0 + c d 0 ) -A kω 0 I m×m -kω 0 I m×m -A M k c x k s x k = B 0 m×p 0 m×p B c u k s u k + c d k s d k , (20)
where I m×m stands for an m × m identity matrix. Equations [START_REF] Halvorsen | Architecture-independent power bound for vibration energy harvesters[END_REF] guarantees that u n , d n and x n satisfy the dynamics (1) at all times, hence representing a consistent dynamic evolution for the system. We introduce the following vectors of Fourier coefficients:

u = [c u 0 , c u 1 , s u 1 . . . , c u n , s u n ] , d = [c d 0 , c d 1 , s d 1 . . . , c d n , s d n ] , x = [c x 0 , c x 1 , s x 1 . . . , c x n , s x n ] , (21)
and express equalities [START_REF] Halvorsen | Architecture-independent power bound for vibration energy harvesters[END_REF] in the following compact form:

M e x = B e (u + d ) , (22)
where M e and B e are sparse block matrices defined as

M e = diag (M 0 , M 1 , . . . , M n ) , B e = diag (B, B, . . . , B) . (23)

Cost function

Expressing the system state, input and excitation as a finite sum of harmonic terms, the cost function of (2) in Problem 1 can be rewritten as a quadratic function of u, x and d which, using [START_REF] Prodanovic | Control and filter design of three-phase inverters for high power quality grid connection[END_REF], can be reduced to u and d only.

In particular, indicating with J the value of the cost function associated to the case of bandlimited dynamics discussed in the previous section (see [START_REF] Prodanovic | Control and filter design of three-phase inverters for high power quality grid connection[END_REF]) and noticing that the harmonic terms in the FS are orthogonal functions in L 2 ([0, T ]) leads to the following equality:

J = T 0 (x n P x n + u n Qu n ) dt = T c x 0 P c x 0 + c u 0 Qc u 0 + + 1 2 n k=1 c x k P c x k + s x k P s x k + c u k Qc u k + s u k Qs u k ( 24
)
Using the coefficient vectors defined in [START_REF] Hals | Constrained optimal control of a heaving buoy wave-energy converter[END_REF], J can be rewritten as follows:

J = T 2 (x D P x + u D Q u) , (25)
where D P and D Q are block-diagonal sparse matrices defined as

D P = diag (2P, P, . . . , P ) , D Q = diag (2Q, Q, . . . , Q) . ( 26
)
Since P and Q are symmetric and positive semidefinite, D P and D Q are also symmetric. Expressing x as a function of u and d through [START_REF] Prodanovic | Control and filter design of three-phase inverters for high power quality grid connection[END_REF], and omitting additive terms and constant positive factors, minimisation of J in ( 25) is equivalent to the minimisation of a quadratic function Ĵ of u, with hessian H and linear term coefficient f defined as follows:

Ĵ(u) = 1 2 u Hu + u f, with H =B e M - e D P M -1 e B e + D Q , f = B e M - e D P M -1 e B e d , (27)
where matrix M e is invertible because A is Hurwitz (a detailed proof of this fact is given in Appendix A).

175

The following important property holds.

Proposition 4

The discretised cost function Ĵ(u) in (27) is convex with respect to u.

Based on (27) we may verify that H is positive semidefinite. Indeed, the first term in the expression of H is positive semidefinite (it is in the form (.) D P (.), with D P positive semidefinite) and the second term D Q is positive semidefinite. Therefore, H is the sum of positive semidefinite matrices, thus

Inequality constraints

A reformulation of the inequality constraints of (2) in Problem 1 in (2) in terms of the problem unknown u is presented in this section.

To this end, we approximate the constraint function g(t) as a finite sum of harmonic terms, in the same fashion as the state, the input, and the excitation:

g(t) ≈ g n (t) = c g 0 + n k=1 c g k cos (kω 0 t) + s g k sin (kω 0 t) , (28)
with c g k ∈ R l , s g k ∈ R l being a set of Fourier coefficients stacked in the following vector:

g = [c g 0 , c g 1 , s g 1 . . . , c g n , s g n ] . (29)
Furthermore, we introduce matrix T n (t), whose elements are functions of time:

T n (t) = I l×l , cos(ω 0 t)I l×l , sin(ω 0 t)I l×l , . . . , cos(nω 0 t)I l×l , sin(nω 0 t)I l×l .

Using ( 29) and (30) and rewriting the constraints of ( 2) in Problem 1 with the finite FS expansions of u and x in [START_REF] Renno | On the optimal energy harvesting from a vibration source[END_REF] we obtain:

T n (t)A n,x x + T n (t)A n,u u + T n (t)g ≤ b, (31)
where A n,x and A n,u are defined as

A n,x = diag (A x , . . . , A x ) , A n,u = diag (A u , . . . , A u ) . (32)
Substituting ( 22) in (31), the constraint can be finally rewritten in terms of u as follows:

A n (t)u + ĝn (t) ≤ b, with

A n (t) = T n (t) A n,x M -1 e B e + A n,u , ĝn (t) = T n (t) A n,x M -1
e B e d + g .

(33)

Frequency-discrete optimisation problem

Based on the recasting of the system state, input, excitation and the cost function discussed above, the LQ optimisation problem in (2) takes the form of a finite-dimensional minimisation problem, as expressed in the following.

Problem 2 Frequency-discrete LQ periodic optimisation problem:

   min u∈R p(2n+1) Ĵ(u) = 1 2 u Hu + u f subj. to: A n (t)u + ĝn (t) ≤ b, ∀t ∈ [0, T ] (34)
Problem 2 is a QP problem in the variable u subject to infinitely many linear parametric constraints to be verified for each value of the parameter, t, within the interval [0, T ]. The solution to Problem 2 represents a suboptimal feasible solution to (2) in Problem 1. In contrast to the solution to Problem 1, the solution to Problem 2 produces a smooth bandlimited evolution of the system, hence providing relevant target optimal dynamics for practical systems, governed by technological constraints.

With this motivation, we hereby focus on the solution of (34) in Problem 2. In principle, under the assumption of H being positive definite, one could reduce the problem into a constrained minimum norm problem over a convex set defined by infinite inequalities. The dual of this problem consists in the maximization of the concave support functional associated to the constraint set over the unit ball, which in principle might seem to be a simpler problem. Besides the strong requirement of strict convexity of the cost function, the computation of that support functional is not trivial, therefore there is no apparent advantage in pursuing this approach.

Problem reformulation with finite-dimensional constraints

Problem 2 is subject to linear continuous-time constraints that need to be verified over a continuous set of time instants t. To enable a computationally efficient solution to the problem, it is hence necessary to further reshape it, e.g., by replacing the infinite dimensional constraints with a finite set of constraints. Regarding this, a numerical approach has been developed in previous works [START_REF] Yen | Linear quadratic optimal control via Fourier-based state parameterization[END_REF][START_REF] Bacelli | Numerical optimal control of wave energy converters[END_REF] that relies on a relaxation of the constraints. According to that approach, the original constraints are only applied in a finite set of time instants, hence introducing a discretisation of the considered time horizon [0, T ] (i.e., one period). The number of the resulting constraints depends on the number of time-samples of the considered time grid, which is independent of the dimension of the unknown variable u. Using a fine discretisation of a time period with a large number of time samples provides a good approximation of the original continuous-time constraints, but it leads to a significant number of constraints (with consequent computational burden) and possible constraint violations in the intersample intervals.

With reference to optimisation (34) in Problem 2, rather than simply using a uniformly spaced grid we adapt a nonuniform grid where we impose a restricted versions of the constraints to get a feasible solution. The corresponding optimisation problem is given next.

Problem 3 Parametric periodic optimisation problem with discretised α-restricted constraints

(P α )    min u∈R p(2n+1) Ĵ(u) = 1 2 u Hu + u f subj. to A n (t k )u + ĝn (t k ) ≤ αb, k = 1, . . . , q (35)
where {t k } q k=0 ⊂ [0, T ] is a finite set of q + 1 time samples, and 0 < α ≤ 1 is a real parameter.

The set T αq = {t k } q k=0 is hereby called the constrained set of time samples (CSTS) for problem (P α ). Since b is positive-valued, reducing α leads to a restriction of the set of the feasible solutions. Nonetheless, if the feasible region for Problem 1 (and, hence, Problem 2) has non empty interior, then there exists a minimum value for α above which the feasible region for (P α ) is non-empty. A special case for the previous problem is α = 1. (P 1 ) has the same objective function as Problem 2 and a wider set of feasible solutions (the inequality constraints are imposed only on a finite set of instants, thus they are less stringent than in the original problem). The optimal solution to (P 1 ) leads to continuous time-series u n and x n that generally violate the constraints of Problem 2 outside of the CSTS (see Fig. 2a).

In contrast, we hereby show that solving (P α ) (with α sufficiently close to 1) potentially allows finding a solution that satisfies to the original finite-dimensional constraints (see Fig. 2b), provided that a suitable CSTS is built.

Some properties of discrete-time constraint approximation

In this section, we discuss the implications of imposing the constraints of Problem 2 over a finite set T αq .

The right-hand side term of the constraint in (33) can be expressed as a finite sum of harmonics 235 as follows:

h(t) = A n (t)u + ĝn (t) = T n (t) A n,x M -1 e B e (d + u) + A n,u u + g , ( 36
)
where the factor into square brackets represents the vector of the Fourier coefficients of h.

The component-wise representation of the constraint in Problem 2 can be thus rewritten in the following form:

h i (t) ≤ b i , with h i (t) = c hi 0 + n k=1 c hi k cos (kω 0 t) + s hi k sin (kω 0 t) , ( 37
)
where h i (t) and b i , with i = 1, . . . , l, are the i-th component of h(t), and b, respectively. With reference to (37), we consider two consecutive time instants t k and t k+1 at which the constraint is respected, and investigate the conditions that prevent constraint violations in the interval [t k , t k+1 ]. Function h i (t) is smooth (it is a finite sum of harmonic terms). Therefore, in case a constraint violation occurs in the interval [t k , t k+1 ], there must exist a point t M in that interval such that ḣi (t M ) = 0 (see Fig. 3), where the violation is maximum. In order to deduce an expression for f i (t M ), we use a Taylor series expansion with Lagrange remainder around t M . Using a first-order Taylor expansion and evaluating it in t k and t k+1 leads to the following expressions:

h i (t k ) = h i (t M ) + ḣi (t ν1 )(t k -t M ), h i (t k+1 ) = h i (t M ) + ḣi (t ν2 )(t k+1 -t M ) ( 38
)
with t ν1 , t ν2 ∈ [t k , t M ] (according to Lagrange theorem). In a similar manner, using a second-order Taylor expansion around t M (where ḣi (t M ) = 0) evaluated at t k and t k+1 leads to the following expressions:

h i (t k ) = h i (t M ) + ḧi (t ν3 ) 2 (t k -t M ) 2 , h i (t k+1 ) = h i (t M ) + ḧi (t ν4 ) 2 (t k+1 -t M ) 2 (39) with t ν3 , t ν4 ∈ [t k , t M ].
Using (37), the following conditions for the first and second order derivatives of h i can be deduced:

| ḣi (t)| ≤ ω 0 n k=1 k (c hi k ) 2 + (s hi k ) 2 = F in , | ḧi (t)| ≤ ω 2 0 n k=1 k 2 (c hi k ) 2 + (s hi k ) 2 = L in . ( 40
)
Using these inequalities, a set of sufficient conditions for (37) to hold over each time interval [t k , t k+1 ] can be found.

Proposition 5 Denoting m ik = min {h i (t k ), h i (t k+1 )} and M ik = max {h i (t k ), h i (t k+1 )}, each one of the following set of inequalities is sufficient for inequality (37) to hold for all t ∈ [t k , t k+1 ]:

m ik + F in • (t k+1 -t k ) ≤ b i , M ik + F in 2 • (t k+1 -t k ) ≤ b i , m ik + L in 2 • (t k+1 -t k ) 2 ≤ b i , M ik + L in 8 • (t k+1 -t k ) 2 ≤ b i . (41)
To prove Proposition 5, we first notice that a sufficient condition ensuring constraint (37) in [t k , t k+1 ] is:

h i (t M ) ≤ b i .

An algorithm for the identification of a suboptimal feasible solution to Problem 2

In this section we present a procedure allowing to iteratively build a feasible solution to (34) in Problem 2, based on a sequence of solutions to different instances of ( 35) in Problem 3. The proposed algorithm relies on a given choice for α < 1 (close to 1) and employs sufficient conditions (41) of Proposition 5 to iteratively build a CSTS such that the chosen solution to (P α ) is feasible for Problem 2.

The steps of the procedure are summarised in Algorithm 1. At each iteration, Algorithm 1 solves (35) in Problem 3 using the current CSTS candidate; sufficient conditions (41) are then tested for each constraint in each interval between consecutive time samples; whenever the test fails, the midpoints of those intervals are added to the CSTS, and the algorithm is iterated until at least a sufficient condition holds in each interval for each of the constraints.

The algorithm clearly provides a solution within the feasible region of Problem 2 in (34), since sufficient conditions (41) guarantee the constraints of Problem 2 over the entire time interval [0, T ]. The outputs of the algorithm are: the CSTS T αq * at the last step, the optimal input coefficients u * α | T αq * , and the optimal value of the cost function J * α | T αq * .

Proposition 6 Let J * be the optimal cost resulting from the solution to (34) in Problem 2, then

lim α→1 J * α | T αq * = J * . (42)
To prove this statement, we consider the following frequency-discrete LQ optimisation problem with α-parametric continuous-time constraints:

   J * α = min u∈R p(2n+1) Ĵ(u) = 1 2 u Hu + u f subect to A n (t)u + ĝn (t) ≤ αb, ∀t ∈ [0, T ] (43)
Eq. ( 43) is a reformulation of (34) in Problem 2 with restricted constraints. The following inequality holds:

J * ≤ J * α | T αq * ≤ J * α , ∀α < 1. ( 44
)
The application of Algorithm 1 leads to a feasible solution to (34) in Problem 2, thereby justifying the left-hand-side inequality. The problem given in (43), in turns, has a smaller feasible region as compared to Problem 2, thus proving the right-hand side inequality in (44). Note also that Problem 2 is a special case of (43) (with α = 1). The problems in (43) are a class Algorithm 1: Computation of a suboptimal feasible solution to (34) in Problem 2 Input: α, T αq 1 ⊇ {0, T } choose α ∈ (0, 1) and an initial CSTS candidate Tαq 1 with 1 + q 1 ≥ 2 samples satisfying to t 0 = 0, tq 1 = T .

j ← 1

Initialise iteration counter repeat u αq j ← solution to (P α ) with T αq j as the CSTS Solve Problem 3

Use u αq j to evaluate h(t) (Eq. ( 36)) at each sample in

T αq j T αq j+1 ← {t 0 } Initialise next-step CSTS for k ← 1 to q j -1 do if ∃ i ∈ {1, . . . ,
q j+1 ← T αq j+1
Update the cardinality of the CSTS

j ← j + 1 until q j = q j-1 return T αq * ← T αq j , u * α | T αq * ← u αq j , J * α | T αq * ← Ĵ u * α | T αq *
of convex problems, whose feasible (convex) set continuously depends on α, and whose optimal cost consequently is a continuous function satisfying

lim α→1 J * α = J * . (45)
Combining this condition with (44) immediately proves Proposition 6.

Proposition 7 Denote by u αq j the solution to (P α ) at the generic j-th step of Algorithm 1. If ||u αq j || 265 is uniformly bounded for all j (namely, ∃ U > 0 : ||u αq j || < U ∀j) then Algorithm 1 converges in a finite number of steps.

Based on (36), since d and g are bounded, if u is also bounded, then the coefficients of h i (t) are bounded. Hence, F in and L in in (40) are also bounded (e.g., ∃ Fin such that F in < Fin at each step). If the algorithm was iterated indefinitely, the lengths of the intervals [t k , t k+1 ] where none of the sufficient conditions (41) hold would decrease and become indefinitely small. In particular, after a finite number of steps, the following condition would hold on each one of the intervals where none of (41) applies:

t k+1 -t k ≤ 2(1 -α)b i Fin , ∀ i = 1, . . . , l. (46)
Noticing that M ik ≤ αb i , (46) immediately implies that the second condition in (41) holds in the above mentioned time intervals, causing the algorithm to stop. Proposition 7 is then proven by contradiction.

• If Ĵ(u) is strictly convex (i.e., H is positive definite), then ||u αq j || is uniformly bounded. The solution to (P α ) at each step of Algorithm 1 provides a sequence of values for the cost function which are upper bounded by the optimal cost for (43):

Ĵ(u αq j ) ≤ J * α . (47)
Being H positive definite, the following inequality holds:

Ĵ(u αq j ) ≥ λ H,min 2 ||u αq j || 2 -||f || • ||u αq j ||, (48)
where λ H,min > 0 is the smallest eigenvalue of H. Combining (47) with (48) immediately proves that ||u αq j || is uniformly bounded.

• It might happen that the constraints in (34) imply that

-ū ≤ u ≤ ū. ( 49
)
In this case the constraints of (35) in Problem 3 may be augmented with (49) to also force the sampled problem to satisfy those uniform bounds for any CSTS candidate, in addition to imposing the hypothesis of Proposition 7. Since the feasible set in Problem 2 is included within the set defined by (49), conditions (42), ( 44) and (45) still hold, in spite of the considered recasting of Problem 3.

Proposition 8 For a given choice of α, denote by T αq * the CSTS resulting from the application of Algorithm 1 to (P α ) in (35), and with J * α | T αq * the resulting optimal cost. Denoting by J * 1 | T αq * the optimal cost for problem (P 1 ) in (35) solved using T αq * as the CSTS, the following inequalities hold:

J * 1 | T αq * ≤ J * ≤ J * α | T αq * , (50)
where J * is the optimal cost for (34) in Problem 2.

The right-hand side inequality in (50) has been previously proven (see (44)). Problem (P 1 ) has the same objective function as (34) in Problem 2, but it has a larger feasible set (i.e., a solution to Problem 2 is always feasible for (P 1 ) but not vice-versa), hence proving the left-hand-side inequality and the proposition.

In practice, solving (P α ) with Algorithm 1 and using the resulting CSTS to solve (P 1 ) provides an upper and a lower bound for the optimal cost of the problem in (34). Eq. (50), combined with Proposition 6, provides a tool to achieving an increasingly accurate estimate of J * , which can be obtained by iteratively applying Algorithm 1 with increasing values of α.

Case study 1 -Mechanical oscillator energy harvester

Consider a direct-force energy harvester (Fig. 1a right) and a harmonic excitation force e = ē cos(ω 0 t), with ω 0 = 2π/T . For the sake of generality, we introduce the following set of dimensionless variables:

t = t T , z = k ē z, ż = kT ē ż, z = kT 2 ē z, u = u ē , e = e ē . (51)
With these definitions, the equations of motion (3) can be rewritten in a dimensionless form:

γ 1 z + γ 2 ż + z = cos(2πt) + u, with γ 1 = M/(kT 2 ), γ 2 = c/(kT ), (52)
and are fully characterised by the dimensionless numbers γ 1 and γ 2 .

Table 1: Numerical parameters of the first case study -mechanical oscillator energy harvester. With reference to (4), we consider a constraint on the maximum dimensionless oscillation amplitude (|z| ≤ z), with no constraints on the control amplitude. With these assumptions, the problem is in the same form as (34) in Problem 2. As an example, we consider the numerical parameters reported in Tab. 1. The considered system has an excitation frequency larger than the natural frequency (i.e., the ratio of the excitation and the natural frequencies, 2πγ 0.5 1 , is greater than 1) and it is underdamped (i.e., the damping ratio 0.5γ -0.5 1 γ 2 is lower than 1). These conditions are representative of some practical applications, e.g., some families of WEC [START_REF] Whittaker | Nearshore oscillating wave surge converters and the development of oyster[END_REF]. Owing to the symmetry of the excitation, the problem can be formulated by omitting zero-order harmonics (i.e., c u 0 = 0) and redefining the matrices of the cost and the constraints (Eq. ( 22) and (33)) by omitting the corresponding zero-order blocks. With this assumption, it can be easily verified that the hessian H of cost (27) (defined as in Appendix B, but without the diagonal block H 0 ) is positive definite, hence guaranteeing convergence of Algorithm 1 in a finite number of steps according to Proposition 7.

We apply Algorithm 1 with α = 0.99 and an initial set of time instants T αq 1 with 5 equally-spaced time samples. The numerical solution of problem (P α ) at the different iterations is tackled using the quadprog command in Matlab. The algorithm converges after 5 iterations. The profiles of the natural variables of the problem (z, e and u) at each iteration are shown in Fig. 4.

It is worth noticing that at the first iterations, due to the large oscillations in the profile of z, the sufficient conditions in (41) do not hold anywhere, and the algorithm uniformly increases the grid of the constrained time samples. In the last few iterations, since the control profile is close to the limiting solution, the algorithm only introduces extra time samples in the intervals where the constraint is active (i.e., |z| z), in order to enforce the sufficient conditions. The number of constrained time samples at convergence is 1 + q * = 36 (see Fig. 4e).

The optimal solution (Fig. 4e) shows a well-known behaviour, referred to as latching [START_REF] Falnes | Ocean waves and oscillating systems: linear interactions including wave-energy extraction[END_REF] in control theory of mechanical energy harvesters with bounded oscillations [START_REF] Halvorsen | Architecture-independent power bound for vibration energy harvesters[END_REF]. This behaviour holds the oscillator fixed at the maximum-displacement positions during certain time intervals, so as to bring the optimal position profile z(t) into phase opposition with respect to the excitation e, mimicking the unconstrained solution.

Based on (50), we used the CSTS T αq * resulting from the application of Algorithm 1 to solve (P 1 ) (i.e., (34) in Problem 2 with α = 1), and we found that the maximum error in the estimate of J * (optimal cost for Problem 2) is J * α | T αq * -J * 1 | T αq * / J * α | T αq * = 0.7%, as proven in Proposition 8. Fig. 4f shows the trend of the amplitudes of the harmonic components of u (i.e., the output of Algorithm 1). Even order terms have an amplitude equal to zero, i.e., the optimal control enjoys a half-wave symmetry, namely u(t + T /2) = -u(t). The amplitude of the odd coefficients decreases with increasing order, exception made for the last term (k = 9) which is larger than the preceding one (k = 7). In practice, the solution algorithm compensates for the truncated number of harmonic components (with respect to the frequency-unbounded solution) through an increase of the highest order term. In Fig. 4e, the optimal solution provided by Algorithm 1 is compared to a numerical solution to (4) 335 obtained using a time domain method with a terminal state constraint [START_REF] Rosati Papini | Control of an oscillating water column wave energy converter based on dielectric elastomer generator[END_REF]. This numerical solution has been obtained using a discretisation of the system dynamics in the time domain (with a fixed time-step ∆t = 5 • 10 -3 ), and using a formulation of the problem similar to the one traditionally used in model predictive control [START_REF] Camacho | Model predictive control[END_REF]. A comparison of the time and the frequency domain solutions shows that, in spite of similar trends and peak values of the controlled force u, the time domain solution has 340 a non-smooth profile with discontinuities in the force values, which are in fact not implementable with practical PTO transducers. The optimal profiles for z, in contrast, are indistinguishable in the two cases: the frequency response of the damped oscillator indeed limits the extent to which higher-order harmonics in the control input can affect the system dynamics.

Case study 2 -Output voltage tracking in DC/AC converters

We consider a single-phase inverter (Fig. 1b) that is required to track a sinusoidal output voltage set-point v r = v cos(ω 0 t). As in the previous example, we define a set of dimensionless variables:

t = t T , q = q vC , q = i = T vC i, q = T 2 vC q, u = u v , v = v v , v r = v r v . (53)
The dynamics [START_REF] Han | A unified numerical scheme for linearquadratic optimal control problems with joint control and state constraints[END_REF] of the circuit can be thus expressed in terms of the following dimensionless form:

γ 1 q + γ 2 q + q = u, with γ 1 = LC T 2 , γ 2 = RC T . (54)
In the following, we investigate the optimal bandlimited input voltage u for the RLC circuit (Eq. ( 13)) subject to an amplitude constraint (u ≤ ū) (resorting to the form (34) in Problem 2). We consider the numerical data reported in Tab. 2, i.e., an underdamped system (0.5γ -0.5 1 γ 2 < 1) that, in contrast with the previous example, has a natural frequency which is larger than the excitation frequency (2πγ 0.5 1 < 1). As in the previous example, we cast the problem as (35) in Problem 3, omitting zero-order harmonic terms and obtaining a positive definite hessian H for cost function (27). We apply Algorithm 1 considering an initial CSTS T αq 1 with 5 equally-spaced time samples. The algorithm converges after 7 iterations. The number of elements in the CSTS T αq * at the last iteration is 1 + q * = 98. The profiles of the dimensionless voltages of the circuit (v, v r and u) at different iterations are shown in Fig. 5. As in the previous example, the algorithm uniformly increases the number of constrained samples during the first iterations (when significant constraint violations occur outside of the constrained set) and then adds samples in the time intervals in which the constraint is active. Since in the considered scenario the condition given by ( 14) is not satisfied, the optimal solution of the problem is non-trivial (i.e., a perfect tracking of the reference trajectory is impossible). Similarly to the previous example, the resulting optimal control has a piecewise-constant profile, assuming the maximum/minimum values (set by the constraint) over finite time intervals.

Following the approximated solution of the problem through Algorithm 1, the resulting CSTS T αq * is used to solve (P 1 ) (see Proposition 8) and to quantify the maximum error in the estimate of J * (optimal cost for (34) in Problem 2), namely: J * α | T αq * -J * 1 | T αq * / J * α | T αq * = 0.2%, as proven in Proposition 8. In Fig. 6a, we report the amplitudes of the harmonic components of u resulting from the application of Algorithm 1. The optimal control input has a half-wave symmetry (i.e., even order harmonics' amplitudes are equal to zero) and the coefficients' amplitude shows a smooth decreasing trend with increasing order.

In Fig. 5d, the discussed solution is compared with that obtained via a time domain approach (using a discretisation of the dynamics with a time-step ∆t = 5 • 10 -3 ). The two numerical procedures (frequency and time domain) provide sensibly different profiles of the input u. In particular, the time domain solution includes some high-frequency oscillations (around t = 0.17 s and t = 0.67 s) between the limit admissible values (|u| = ū), which do not appear in the frequency domain solution as a consequence of the finite number of harmonics. Such differences in the input, however, do not generate significant differences in terms of the system output v and its distance from the reference trajectory. Solving the frequency domain problem (via Algorithm 1) with a larger number of harmonics leads to a solution with high-frequency oscillations in the optimal profile of u, similarly to the optimal time domain solution. As an example, Fig. 6b compares the numerical solutions of the problem in the case of a time domain approach and a frequency domain approach with maximum order of the Fourier coefficients n = 60 (as opposed to n = 10 used in Fig. 5). This example clearly shows that the frequency domain approach provides a smoother solution for the control input (easier to pursue in practice) as compared to the time domain approach, while preserving the features of the time domain optimal solution in terms of the system output response.

Conclusions

We addressed dynamic programming in periodic systems via a functional approach based on the Fourier series (FS). As compared to time domain methods, which rely on a discretisation of the dynamics through the state and input evaluation over a finite grid of time instants, frequency domain methods based on the FS provide smooth optimal profiles that can be practically employed to develop control heuristics for real systems. Motivated by two technical applications, we addressed linear quadratic problems on periodic steadystate responses. We provided a general formulation that allows casting the optimisation problem as an equivalent minimisation problem with parametric constraints, whose unknowns are the coefficients of a FS with a finite number of elements. Solving this minimisation problem leads to smooth bandlimited control inputs. As compared to previous works, this paper extends the framework of Fourier-based optimisation methods in a rigorous way by carefully discussing the implications of a convex frequency domain formulation for the periodic problem, and proposes a new numerical procedure providing approximated globally-feasible solutions. The proposed procedure relies on the iterative solution of a set of quadratic problems, whose constraints are built by sequentially restricting the constraints of the original optimisation problem and relaxing the set of time instants in which such constraints are enforced. Two case studies have been discussed, which show that the proposed framework can be applied to solve problems of actual interest in research and engineering (optimal control of vibrating energy harvesters and of DC/AC converters) which can be represented in a simplified and essential way as periodic LQ problems.

Figure 1 :

 1 Figure 1: Application examples. (a) Lumped parameter models of oscillating energy harvesters: inertial generator (left), and direct-force generator (right). (b) Circuit schematization of a single-phase DC/AC converter.

  145

Figure 2 :

 2 Figure 2: (a) Approximating Problem 2 with (P 1 ) leads to constraint violation. (b) Solving (P α ) provides a feasible solution to Problem 2 for an appropriate choice of the CSTS.

  240

Figure 3 :

 3 Figure 3: Constraint violation in a time interval between two points t k and t k+1 where constraint (37) is respected.

Figure 4 :

 4 Figure 4: (a-e) Profiles of z, e and u at each iteration (1 to 5) of Algorithm 1. The green dashed lines represent the constraint level z ≥ |z|. The asterisks on the time axis represent the elements of the CSTS at each iteration. Figure e) also reports a non-smooth numerical solution obtained with a time domain formulation with no restrictions on the control bandwidth (purple dashed line).Figure f) shows the amplitude of the ten harmonic coefficients in the FS of the output u of Algorithm 1.

Figure 5 :

 5 Figure 5: (a-d) Profiles of z, e and u at different iterations (1 to 7) of Algorithm 1. The green dashed lines represent the constraint level ū ≥ |u|. The asterisks on the time axis represent the elements of the CSTS at the different iterations. Figure d) also reports a numerical solution to the problem obtained from a time domain formulation with no restrictions on the control bandwidth (purple dashed line).

Figure 6 :

 6 Figure 6: (a) Amplitude of the harmonic terms in the FS of u (output of Algorithm 1). (b) Optimal profiles of z, e and u obtained using Algorithm 1 with n = 60 harmonics.

  Following this, the first and the third conditions in (41) are easily obtained from Eqs. (38-39) replacing the derivatives with their upper-bounds in (40) and noticing that (t M -t k ), (t k+1 -t M ) ≤ (t k+1 -t k ). The second and the fourth conditions in (41) are obtained in a similar manner, noticing that min{t M -t k , t k+1 -t M } ≤ (t k+1 -t k )/2.

  l} s.t. none of suff. cond. (41) holds then T αq j+1 ← T αq j+1 ∪ t k + t k+1 2 , t k+1 Include interval midpoint in the CSTS else T αq j+1 ← T αq j+1 ∪ {t k+1 }

	Do not include the midpoint in the CSTS
	end
	end

Table 2 :

 2 Numerical parameters of the second case study -voltage tracking in inverters.

	γ 1	γ 2	ū	n	α
	8 • 10 -3 0.06 0.55 10 0.99

immediately proving Proposition 4.

Whether or not u satisfies the boundedness assumption in Proposition 7 depends on the problemspecific form of the cost function and the constraints. Two scenarios in which the boundedness condition in Proposition 7 holds are as follows.

Case studiesIn this section, we present two numerical case studies involving the applications discussed in Sect. 2.1 (mechanical energy harvester) and Sect. 2.2 (inverter control). With reference to a sinusoidal excitation d, bandlimited solutions are investigated via the application of Algorithm 1 and the sandwich bound in (50).

Acknowledgement

GM received funds from Regione Toscana (Italy) via project EOLO (FAR FAS 2014-A). LZ received funds from the Agence Nationale de la Recherche (ANR) via grant "Hybrid And Networked Dynamical sYstems" (HANDY), number ANR-18-CE40-0010.

APPENDICES A Proof of the invertibility of M e

We hereby prove that if A is a Hurwitz matrix (i.e., its eigenvalues have strictly negative real part) then M e is invertible. M e is a block-diagonal matrix defined as in [START_REF] Mariéthoz | Explicit model-predictive control of a PWM inverter with an LCL filter[END_REF]. Since the eigenvalues of A are different from zero, M 0 is invertible (see [START_REF] Halvorsen | Architecture-independent power bound for vibration energy harvesters[END_REF]). Matrices M k are invertible provided that the following condition holds for all k = 1, . . . , n:

as this ensures the existence of the inverse matrices for M k :

Since A is Hurwitz, the eigenvalues of A 2 (i.e., the squares of the eigenvalues of A) are either positive real numbers or complex numbers, hence automatically satisfying condition (A.1).

Finally, M e is invertible since each of its diagonal blocks is invertible.

B Explicit expression for H

Based on (27), matrix H is block diagonal and it has the explicit expression given by (B.1), where the first diagonal block (H 0 ) belongs to R p×p , and the other blocks (H k ) belong to R 2p×2p .